##Analysis file for Surrogacy experiment##
library(tidyverse)
library(ggplot2)
library(gganimate)
library(readr)
library(broom)
library(estimatr)
library(interactions)
library(jtools)
library(ggplot2)
library(ggstance)
library(rdrobust)   
library(rddensity) 
library(patchwork)
library(ggridges)
library(colorspace)
library(ggpubr)
library(modelsummary)
library(interplot)
library(margins)
library(ggdist)

surrogacy <- read_csv("surrogacyDATA.csv")

surrogacy <- surrogacy %>% 
  mutate(treatment = as.factor(treatment),
         homo = as.factor(homo),
         homo_num = as.numeric(homo),
         prestige = as.factor(prestige),
         gender1 = as.factor(gender_nolab),
         education = as.factor(education),
         degree1 = as.factor(degree_nolab),
         race1 = as.factor(race_nolab),
         income1 = as.numeric(income_nolab),
         sexuality1 = as.factor(sexuality_nolab),
         sexactive1 = as.factor(sexactive_nol),
         ideology1 = as.factor(ideology_nol),
         vote20191 = as.factor(vote19_nol),
         voteleave1 = as.factor(voteleave_nol),
         employ1 = as.factor(employ_nol),
         civilstat1 = as.factor(civilstat_nol),
         religion1 = as.factor(relig_nol),
         property1 = as.factor(property_nol))


surr_nocontrol<-surrogacy %>% 
  filter((treatment!=0))





###FIGURE 1###

col1 <- c("#ffb3ba", "#ffdfba", "#ffffba", "#baffc9", "#bae1ff")

model1 <-lm (support ~ treatment, data=surrogacy)
surrogacy$predicted<-predict(model1, surrogacy)

effect_plot(model = model1, pred = treatment, robust=TRUE,
            cat.geom="point", cat.interval.geom="linerange",
            colors="slategrey", cat.pred.point.size=3)+
  labs(title = "Mean level of support by treatment condition",
       subtitle = "Pr(Supports surrogacy)",
       caption= "Treatment group outcomes statistically distinct at p<0.05(**) and p<0.01(***)")+
  ylab("")+
  xlab("Treatment condition")+
  geom_jitter(data=surrogacy, aes(x=treatment, y=predicted, fill=treatment),
              height=.1, width=.45,size=3, alpha=.3,
              pch=21, color="grey57")+
  scale_color_manual(values = col1) +
  scale_fill_manual(values = col1)+ 
  scale_x_discrete(labels=c("0" = "Control", "1" = "Celebrity\ngay couple", "2" = "Celebtrity\nstraight couple",
                            "3" = "Ordinary\ngay couple", "4" = "Ordinary\nstraight couple"))+
  theme(legend.position = "none",
        axis.text.x = element_text(face=2))+
  geom_bracket(
    xmin = c("0", "1"), xmax = c("3", "3"),
    y.position = c(.61, .64), label = c("***", "***"),
    tip.length =-0.02,
    color="slategrey")+
  geom_bracket(
    xmin = c("2", "3"), xmax = c("3", "4"),
    y.position = c(.85, .9), label = c("**", "***"),
    tip.length =0.02,
    color="slategrey")
ggsave("Figure1.tiff",height=16, width=18, units="cm", dpi=400)



##FIGURE 2##
modelinteractX<-lm (support ~ homo*prestige, data=surr_nocontrol)
modelmargX <-
  modelinteractX %>%
  margins(at = list(prestige = c("0", "1"))) %>%
  summary %>%
  as.data.frame() %>%
  filter(factor == "homo1")

modelmargX2 <-
  modelinteractX %>%
  margins(at = list(homo = c("0", "1"))) %>%
  summary %>%
  as.data.frame() %>%
  filter(factor == "prestige1")

ameIplot1<- ggplot(modelmargX, aes(prestige, AME)) +
  geom_point(colour="#C5407EFF", size=2) +
  geom_errorbar(aes(ymax = lower, ymin = upper), size=1, width = 0, colour="#C5407EFF") +
  geom_hline(yintercept = 0, linetype = "dashed", colour="#205C8A") +
  xlab("") +
  ylab("") +
  ylim(-.2,.2)+
  annotate(
    geom="text", x =1.25, y = -.13, size = 3, color = "#C5407EFF", fontface=2,
    label = "-13\n%-points")+
  annotate(
    geom="text", x =1.75, y = 0.03, size = 3, color = "#C5407EFF", fontface=2,
    label = "0.8\n%-points")+
  scale_x_discrete(labels=c("0" = "Ordinary\ncouple", "1" = "Celebrity\ncouple"))+
  labs(subtitle="Conditional AME of (homo)sexuality")+
  theme_minimal()+
  theme(plot.subtitle = element_text(hjust = 0, face="bold"))

ameIplot2<- ggplot(modelmargX2, aes(homo, AME)) +
  geom_point(colour="#C5407EFF", size=2) +
  geom_errorbar(aes(ymax = lower, ymin = upper), size=1, width = 0, colour="#C5407EFF") +
  geom_hline(yintercept = 0, linetype = "dashed", colour="#205C8A") +
  xlab("") +
  ylab("") +
  ylim(-.2,.2)+
  annotate(
    geom="text", x =1.25, y = -.07, size = 3, color = "#C5407EFF", fontface=2,
    label = "-6.7\n%-points")+
  annotate(
    geom="text", x =1.75, y = .07, size = 3, color = "#C5407EFF", fontface=2,
    label = "6.9\n%-points")+
  scale_x_discrete(labels=c("0" = "Straight\ncouple", "1" = "Homosexual\ncouple"))+
  labs(subtitle="Conditional AME of parasocial contact")+
  theme_minimal()+
  theme(plot.subtitle = element_text(hjust = 0, face="bold"))

plotM1<- ameIplot1+ameIplot2+
  plot_annotation(title = "Conditional effect of treatment factors",
                  theme = theme(plot.title = element_text(size = 16, face="bold")))

ggsave("Figure 2.tiff", dpi=300, height=16, width=18, units="cm")



##FIGURE 3##
library(viridisLite)
colors <- viridis(option = "plasma", 
                  begin = 0.15, 
                  end = 0.8, 
                  direction = -1, 
                  n = 5)

model3 <-lm (support ~ treatment*age, data=surrogacy)
surrogacy$predictINT<-predict(model3, surrogacy)

control <- subset(surrogacy, treatment== 0)
treat1 <- subset(surrogacy, treatment== 1)
treat2 <- subset(surrogacy, treatment== 2)
treat3 <- subset(surrogacy, treatment== 3)
treat4 <- subset(surrogacy, treatment== 4)

ggplot(surrogacy, aes(age, support, shape = treatment, group = treatment, color=treatment)) +
  scale_color_manual(values = colors, labels=("0" = "Control")) +
  scale_fill_manual(values = colors) +
  stat_smooth(method = "lm_robust", fullrange = TRUE, se=FALSE) +
  theme_minimal()+
  xlim(18,100)+
  labs(y="Pr(Support surrogacy)", x="Age",
       title="Conditional ATE across age distribution")+
  geom_jitter(data=control, aes(x=age, y=predictINT, fill=treatment), 
              alpha=.2, height=.03, width=2,
              pch=21, color="#FCA636FF")+
  geom_jitter(data=treat1, aes(x=age, y=predictINT, fill=treatment), 
              alpha=.4, height=.03, width=2,
              pch=21, color="#E77059FF")+
  geom_jitter(data=treat2, aes(x=age, y=predictINT, fill=treatment), 
              alpha=.2, height=.03, width=2,
              pch=21, color="#C5407EFF")+
  geom_jitter(data=treat3, aes(x=age, y=predictINT, fill=treatment), 
              alpha=.2, height=.03, width=2,
              pch=21, color="#9511A1FF")+
  geom_jitter(data=treat4, aes(x=age, y=predictINT, fill=treatment), 
              alpha=.2, height=.03, width=2,
              pch=21, color="#5601A4FF")+
  theme(legend.position = "none",
        axis.text.x = element_text(face=2, color="black"),
        axis.ticks.y = element_blank(),
        axis.text.y = element_text(face=2, color="black"),
        axis.title.x = element_text(face=2),
        axis.title.y = element_text(face=2),
        plot.title = element_text(face=2))+
  annotate(
    geom="text", x = 85, y =.75, size = 3, color = "#FCA636FF", fontface=2, hjust = 0,
    label = "Slope for\ncontrol group")+
  annotate(
    geom="text", x = 70, y =.85, size = 3, color = "#5601A4FF", fontface=2, hjust = 0,
    label = "Slope for hetero.\ncivilian treatment")+
  annotate(
    geom="text", x = 25, y =.6, size = 3, color = "#9511A1FF", fontface=2, hjust = 0,
    label = "Slope for homo.\ncivilian treatment")+
  annotate(
    geom="text", x = 50, y =.45, size = 3, color = "#C5407EFF", fontface=2, hjust = 0,
    label = "Slope for hetero.\ncelebrity treatment")+
  annotate(
    geom="text", x = 60, y =.3, size = 3, color = "#E77059FF", fontface=2, hjust = 0,
    label = "Slope for homo.\ncelebrity treatment")+
  annotate(
    geom = "curve", x = 90, y = .72, xend = 85, yend = .61, 
    curvature = -.2, arrow = arrow(length = unit(2, "mm")), colour="#FCA636FF")+
  annotate(
    geom = "curve", x = 71, y = .85, xend = 65, yend = .79, 
    curvature = .2, arrow = arrow(length = unit(2, "mm")), colour="#5601A4FF")+
  annotate(
    geom = "curve", x = 35, y = .63, xend = 45, yend = .7, 
    curvature = -.2, arrow = arrow(length = unit(2, "mm")), colour="#9511A1FF")+
  annotate(
    geom = "curve", x = 67, y = .45, xend = 77, yend = .52, 
    curvature = .3, arrow = arrow(length = unit(2, "mm")), colour="#C5407EFF")+
  annotate(
    geom = "curve", x = 76, y = .30, xend = 85, yend = .46, 
    curvature = .2, arrow = arrow(length = unit(2, "mm")), colour="#E77059FF")
ggsave("Figure3.tiff", dpi=400)



###GENDER###
men<-surrogacy %>% 
  filter((gender_nolab==0))

women<-surrogacy %>% 
  filter((gender_nolab==1))

modelX <-lm (support ~ treatment, data=men)
men$predictedX<-predict(modelX, men)
summ(modelX)
modelX2 <-lm (support ~ treatment, data=women)
women$predictedX2<-predict(modelX2, women)
summ(modelX2)


models <- list()
models [['Men']] <- lm (support ~ treatment, data=men)
models [['Women']] <- lm (support ~ treatment, data=women)
modelsummary(models, star=TRUE, output="latex")

colM <- c("#00695C", "#00796B", "#00897B", "#009688", "#26A69A")
colF <- c("#AD1457", "#C2185B", "#D81B60", "#E91E63", "#EC407A")


X1<- effect_plot(model = modelX, pred = treatment, robust=TRUE,
                 cat.geom="point", cat.interval.geom="linerange",
                 colors="slategrey", cat.pred.point.size=3)+
  theme_minimal()+
  labs(
    subtitle = "Men")+
  ylab("")+
  ylim(.5, 1)+
  xlab("")+
  geom_jitter(data=men, aes(x=treatment, y=predictedX, fill=treatment),
              height=.1, width=.4,size=2, alpha=.35,
              pch=21, color="grey57")+
  scale_color_manual(values = col1) +
  scale_fill_manual(values = col1)+ 
  scale_x_discrete(labels=c("0" = "Control", "1" = "Gay\ncelebrity", "2" = "Stright\ncelebtrity",
                            "3" = "Gay\ncivilian", "4" = "Straight\ncivilian"))+
  theme(legend.position = "none",
        axis.text.x = element_text(face=2, color="black"),
        axis.text.y = element_text(face=2, color="black"),
        plot.subtitle = element_text(face="bold", color="black"))+
  geom_bracket(
    xmin = c("0", "3"), xmax = c("1", "4"),
    y.position = c(.95, .92), label = c("**", "***"),
    tip.length =0.02,
    color="black")+
  geom_bracket(
    xmin = c("0", "0"), xmax = c("2", "3"),
    y.position = c(.61, .54), label = c("***", "***"),
    tip.length =-0.02,
    color="black")+  
  annotate(
    geom="text", x = 5, y = .51, size = 3,
    color = "black",
    label = "N=923")

X2<- effect_plot(model = modelX2, pred = treatment, robust=TRUE,
                 cat.geom="point", cat.interval.geom="linerange",
                 colors="slategrey", cat.pred.point.size=3)+
  theme_minimal()+
  labs(
    subtitle = "Women")+
  ylab("")+
  ylim(.5, 1)+
  xlab("")+
  geom_jitter(data=women, aes(x=treatment, y=predictedX2, fill=treatment),
              height=.1, width=.4,size=2, alpha=.35,
              pch=21, color="grey57")+
  scale_color_manual(values = col1) +
  scale_fill_manual(values = col1)+ 
  scale_x_discrete(labels=c("0" = "Control", "1" = "Gay\ncelebrity", "2" = "Stright\ncelebtrity",
                            "3" = "Gay\ncivilian", "4" = "Straight\ncivilian"))+
  theme(legend.position = "none",
        plot.subtitle = element_text(face="bold", color="black"),
        axis.text.y = element_blank(),
        axis.text.x = element_text(face=2, color="black"),
        axis.ticks.y = element_blank())+
  geom_bracket(
    xmin = c("0"), xmax = c("4"),
    y.position = c(.97), label = c("*"),
    tip.length =0.02,
    color="black")+
  geom_bracket(
    xmin = c("3"), xmax = c("4"),
    y.position = c(.65), label = c("***"),
    tip.length =-0.02,
    color="black")+  
  annotate(
    geom="text", x = 5, y = .51, size = 3,
    color = "black",
    label = "N=946")

X1+X2+ 
  plot_annotation(title = 'Conditional ATE across gender subsamples',
                  caption="Treatment group outcomes statistically distinct at p<.1(*), p<0.05(**), & p<0.01(***)",
                  theme = theme(plot.title = element_text(size = 14, face="bold")))
ggsave("Figure 4.tiff",height=16, width=18, units="cm", dpi=400)


